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(57) A method is disclosed for determining the hor- 
izontal resistivity, vertical resistivity and dip of forma- 
tions surrounding a bore-hole, comprising the following 
steps: (a) suspending a logging device in the borehole; 
(b) transmitting electromagnetic energy from a transmit- 
ter location on the logging device, and receiving the 
transmitted electromagnetic energy at receiver loca- 
tions on the logging device for a first transmitter-to-re- 
ceivers spacing associated with the transmitter and re- 
ceiver locations; (c) determining, from the received elec- 
tromagnetic energy, measurement characteristics asso- 
ciated with the first transmitter-to-receivers spacing; (d) 
repeating steps (b) and (c) for a plurality of further trans- 
mitter-to-receivers spacings to obtain measurement 
characteristics for the plurality of further transmitter-to- 
receivers spacings; (e) generating a model of aniso- 
tropic for-mations that has horizontal resistivity Rh, ver- 
tical resistivity Rv, and a Dip angle with respect to a bore- 
hole reference; (f) selecting initial model values of Rh, 
Rv and Dip; (g) computing an error value from the dif- 
ferences, for each of a plurality of the transmitter-to-re- 
ceivers spacings, between the measurement character- 
istics and model measurement characteristics obtained 
from the model using the model values of Rh, Rv, and 
Dip; (h) modifying the model values of Rh, Rv and Dip; 
(i) repeating steps (g) and (h) to minimize the error val- 
ue; and (j) outputting the ultimately modified model val- 
ues of Rh ; Rv and Dip. 
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Description 

This invention relates to the field of earth borehole logging and, more particularly, to a method and apparatus for 
determining properties of anisotropic earth formations. Properties that can be determined using the invention include 
s horizontal conductivity, vertical conductivity, horizontal dielectric constant, vertical dielectric constant, and dip angle. 

Field of the Invention 

The well logging art has recognized for many years that formations surrounding an earth borehole can be aniso- 

10 tropic with regard to conduction of electrical currents (see e.g. K.S. Kunz et al., "Some Effects Of Formation Anisotropy 
On Resistivity Measurements In Boreholes", Geophysics , Vol. 23, No. 4, 1958). The phenomenon of electrical anisot- 
ropy is generally visualized in one of two ways, or a combination thereof, as follows. 

In many sedimentary strata, electric current flows more easily in a direction parallel to the bedding planes than 
transversely to them. A reason for this anisotropy is that a great number of mineral crystals possess a flat or elongated 

15 shape (e.g. mica and kaolin). At the time they were laid down, they naturally took an orientation parallel to the plane 
of sedimentation. The interstices in the formations are, therefore, generally parallel to the bedding plane, and the 
current is able to travel with facility along these interstices which often contain electrically conductive mineralized water. 
Such electrical anisotropy, sometimes called microscopic anisotropy, is observed mostly in shales. 

If a cylindrical sample is cut from a formation, parallel to the bedding planes, the resistivity of this sample measured 

20 with a current flowing along its axis is called the longitudinal (or horizontal) resistivity R h . The inverse of R h is the 
horizontal conductivity, a h . If a similar cylinder is cut perpendicular to the bedding planes, the resistivity measured with 
a current flowing along its axis is called the transversal (or vertical) resistivity F^. The inverse of R v is the vertical 
conductivity, a v . The anisotropy coefficients, by definition, is equal \ojR ] /R h (or Jo^/a^ ) . Laboratory measurements 
have shown that X may range from 1 to about 2.5 in different shales. 

25 Furthermore, the formations are often made up of a series of relatively thin beds having different lithologic char- 

acteristics and, therefore, different resistivities (as, for example, sequences of thin shales and hard streaks). In well 
logging systems the distances between the electrodes or antennas are great enough that the volume involved in a 
measurement may include several such thin beds. Since, in this situation, the current flows more easily along the more 
conductive streaks than transversely to the series of beds, there is effective anisotropy. The effects on resistivity meas- 

30 urements of this "macroscopic" anisotropy are cumulative with the effects of the anisotropy due to the above-described 
microscopic structure of the sediments. Reference can also be made to J.H. Moran et al., "Effects Of Formation Ani- 
sotropy On Resistivity Logging Measurements", Geophysics , Vol. 44, No. 7, 1979, and to R. Chemali etal., "The Effect 
Of Shale Anisotropy On Focused Resistivity Devices", SPWLA Twenty-Eighth Annual Logging Symposium, 1987. 
The determination of R v as well as R h can be useful in various situations. For example, consider the case where 

35 the formation consists substantially of two types of material with resistivities R-| and R 2 with respective volume fractions 
a and 1-a. The effective horizontal and vertical resistivities R h 

R 1 R 2 

4Q h (a R 2 ) + ((1 - ajRj) V } 

and R v are given by 



45 



R v =(aR 1 )^((1-a)R 2 ) (2) 



If a is known, such as in a shale sequence where a gamma ray measurement or a spontaneous potential measurement 
has been used to provide the shale fraction, R-, and R 2 can be determined from (1) and (2) if R h and R v are known. 
In situations where the borehole intersects the formations substantially perpendicular to the bedding planes, con- 

50 ventional induction and propagation well logging tools are sensitive almost exclusively to the horizontal components 
of the formation resistivity. When the borehole intersects the bedding planes at an angle, the tool readings contain an 
influence from the vertical resistivity as well as the horizontal. This is particularly true when the angle between borehole 
and the normal to the bedding planes becomes large, such as in directional or horizontal drilling where angles near 
90° are commonly encountered. In these cases, the influence of vertical resistivity can cause discrepancies between 

55 measurements in these wells and measurements taken of the same formation in nearby vertical wells, thereby pre- 
venting useful comparison of these measurements. In addition, since reservoir evaluation is typically based upon data 
from vertical wells, use of data from wells drilled at high angles may produce erroneous estimates of formation pro- 
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ducibility if proper account is not taken of the anisotropy effect. 

A number of techniques have been proposed for measuring formation anisotropy and/or vertical conductivity, such 
as by providing transmitter and/or receiver coils that are perpendicular to the borehole axis in addition to coils having 
conventional orientations. Reference can be made, for example, to U.S. Patents 4 ; 302,722, 4,302,723, and 4,980,643. 

s Equipment and techniques that determine horizontal and vertical conductivity (or anisotropy) by employing special 

equipment dedicated specifically to such purpose result in increased equipment cost and increased logging time and/ 
or cost. U.S. Patent 5,329,448 discloses a technique for determining horizontal and vertical conductivity (or anisotropy 
determinable therefrom) using measurements that are often available from conventional types of equipment utilized 
for logging earth boreholes during drilling or by wireline. In the '448 Patent, first and second formation conductivity 

10 values are derived from measurements taken in a borehole. The derived first and second formation conductivity values 
are obtained from measurements which are affected differently by the vertical and horizontal conductivities of the 
formations. These may be, for example, conductivities respectively obtained from phase and attenuation measure- 
ments from a receiver pair of a propagation logging device such as that disclosed in U.S. Patent 4,899,112, or con- 
ductivities respectively obtained from the resistive and reactive components of an induction logging apparatus. Hori- 

15 zontal and vertical conductivity model values are selected, and error values are computed from the differences between: 
(i) first and second composite conductivity values computed as a function of the horizontal and vertical conductivity 
model values, and (ii) the first and second derived formation conductivity values. In iterative fashion, the model values 
are modified in a manner that tends to reduce the error values. The '448 Patent continues the process until a prede- 
termined criterion of the error is met, whereupon the modified horizontal and vertical conductivity model values can be 

20 read out as the determined values. 

It is among the objects of the present invention to provide a method and apparatus for determining properties of 
anisotropic earth formations that improves on prior techniques for determination of horizontal and vertical conductivity 
and is also capable of providing dip angle and anisotropic dielectric constant characteristics. 

25 SUMMARY OF THE INVENTION 

In an embodiment of the invention, the dip angle, as well as the horizontal and vertical resistivities, can be deter- 
mined, and it is not necessary to have a priori knowledge of the dip angle in order to determine horizontal resistivity 
and vertical resistivity. In the present application, unless otherwise specified, any references to determination or use 

30 of resistivity are intended to generically mean conductivity as well, and vice versa. These quantities are reciprocals, 
and mention of one or the other herein is for convenience of description, and not intended in a limiting sense. 

In accordance with an embodiment of the invention there is disclosed a method for determining the horizontal 
resistivity, vertical resistivity, and dip of formations surrounding a borehole, comprising the following steps: (a) sus- 
pending a logging device in the borehole; (b) transmitting electromagnetic energy from a transmitter location on the 

35 logging device, and receiving the transmitted electromagnetic energy at receiver locations on the logging device for a 
first transmitter-to-receivers spacing associated with the transmitter and receiver locations; (c) determining, from the 
received electromagnetic energy, measurement characteristics associated with the first transmitter-to-receivers spac- 
ing; (d) repeating steps (b) and (c) for a plurality of further transmitter-to-receivers spacings to obtain measurement 
characteristics for said plurality of further transmitter-to-receivers spacings; (e) generating a model of anisotropic for- 

40 mations that has horizontal resistivity Rh, vertical resistivity Rv, and a Dip angle with respect to a borehole reference; 
(f) selecting initial model values of Rh, Rv and Dip; (g) computing an error value from the differences, for each of a 
plurality of said transmitter-to-receivers spacings, between said measurement characteristics and model measurement 
characteristics obtained from the model using the model values of Rh, Rv, and Dip; (h) modifying the model values of 
Rh, Rv and Dip; (i) repeating steps (g) and (h) to minimize the error value; and (j) outputting the ultimately modified 

45 model values of Rh, Rv, and Dip. 

In a disclosed embodiment, the step of selecting initial model values of Rh, Rv, and Dip comprises: producing 
discretized model values of Rh, Rv and Dip; comparing the measurement characteristics with model measurement 
characteristics obtained using a multiplicity of combinations of discretized model values; and selecting, as the initial 
model values, discretized model values based on the comparisons of the comparing step. This search for the initial 

50 model values is global in nature and is used to provide candidate model values that are close enough to the correct 
final solution so that the next more refined process of inversion is likely to converge to the correct final solution rather 
than to another solution made possible by the complex (and not always orderly) nature of the relationships between 
the measured resistivities at various spacings and the anisotropic formation dip angle. 

Further features and advantages of the invention will become more readily apparent from the following detailed 

55 description when taken in conjunction with the accompanying drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 is a diagram, partially in blockform, of a system in which an embodiment of the invention can be employed, 
and which can be used to practice an embodiment of the method of the invention. 
s Figure 2 is a diagram, partially in block form, of electronics in the measuring apparatus and the acquisition and 

processor electronics of the Figure 1 embodiment. 

Figure 3 is a diagram illustrating a transverse isotropic model with a dip angle. 

Figure 4 is a graph of phase-shift, phase-average, and attenuation resistivity responses as a function of relative 
dip angle for various transmitter-to-receivers spacings, TR, in a transverse isotropic medium with a relatively high 
10 anisotropy ratio. 

Figure 5 is a graph of phase-shift, phase-average, and attenuation resistivity responses as a function of relative 
dip angle for various transmitter-to-receivers spacings ; TR, in a transverse isotropic medium with a relatively low ani- 
sotropy ratio. 

Figure 6, which includes Figure 6B placed below Figure 6A, is a flow diagram which, in conjunction with the flow 
15 diagrams referred to therein, can be used in programming a processor or processors in implementing an embodiment 
of the invention. 

Figure 7, which includes Figure 7B placed below Figure 7A, is a flowdiagram of a routine for controlling a processor 
to implement borehole logging and storage of measurements. 

Figure 8, which includes Figure 8C placed below Figure 8B placed below Figure 8A, is a flowdiagram of a routine 
20 for implementing an initial search to obtain candidates as initial model values. 

Figure 9, which includes Figure 9B placed below Figure 9A, a flow diagram of a routine for implementing inversion 
to obtain model parameter values of Rh, Rv, and Dip. 

Figure 10 is a plot of fractional error, for approximated vertical conductivity and dielectric constant, at various dip 
angles. 

25 

DESCRIPTION OF THE PREFERRED EMBODIMENTS 

Referring to Figure 1. there is illustrated an embodiment of the invention in the form of a logging-while-drilling 
apparatus and method. A platform and derrick 1 0 are positioned over a borehole 1 1 that is formed in the earth by rotary 

30 drilling. A drill string 12 is suspended within the borehole 11 and includes a drill bit 15 at its lower end. The drill string 
12, and the drill bit 1 5 attached thereto, is rotated by a rotating table 1 6 (energized by means not shown) which engages 
a kelly 17 at the upper end of the drill string 12. The drill string 12 is suspended from a hook 18 attached to a traveling 
block (not shown). The kelly 17 is connected to the hook 18 through a rotary swivel 1 9 which permits rotation of the 
drill string 12 relative to the hook 18. Drilling fluid or mud 26 is contained in a pit 27 in the earth. A pump 29 pumps the 

35 drilling fluid 26 into the drill string 12 via a port in the swivel 19 to flow downward through the center of drill string 12. 
The drilling fluid 26 exits the drill string 12 via ports in the drill bit 15 and then circulates upward in the region between 
the outside of the drill string 12 and the periphery of the borehole 11. As is well known, the drilling fluid 26 thereby 
carries formation cuttings to the surface of the earth, and the drilling fluid 26 is returned to the pit 27 for recirculation. 
The small arrows in Figure 1 illustrate the typical direction of flow of the drilling fluid 26. 

40 Mounted within the drill string 12, preferably near the drill bit 15, is a downhole sensing, processing, storing and 

transmitting subsystem 100. Subsystem 100 includes a measuring apparatus 200 of the general type disclosed in U. 
S. Patent 5,594,343 which, in the present embodiment, includes five transmitting antennas T1 , T2, T3, T4 and T5 and 
receiving antennas R1 and R2, and operates in the manner to be described below. The antennas can be of the type 
described in U.S. Patent 4,899,112; that is, coils wound on mounting material around a special section of metal drill 

45 collar which comprises part of the subsystem 100. A communications transmitting portion of the downhole subsystem 
100 includes an acoustic transmitter 56, which generates an acoustic signal in the drilling fluid 26 that is representative 
of the measured downhole conditions. One suitable type of acoustic transmitter, which is known in the art, employs a 
device known as a "mud siren" which includes a slotted stator and a slotted rotor that rotates and repeatedly interrupts 
the flow of drilling fluid to establish a desired acoustic wave signal in the drilling fluid. The generated acoustic mud 

50 wave travels upward in the drilling fluid through the center of the drill string at the speed of sound in the drilling fluid. 
The acoustic wave is received at the surface of the earth by transducers represented by reference numeral 31 . The 
transducers, which are, for example, piezoelectric transducers, convert the received acoustic signals to electronic 
signals. The output of the transducers 31 is coupled to the uphole receiver subsystem 90 which is operative to demod- 
ulate the transmitted signals, which are then coupled to processor 85 and recorder 45. 

55 Transmitter 56 is controlled by transmitter control and driving electronics 57 which includes analog-to-digital (A/ 

D) circuitry that converts the signals representative of downhole conditions into digital form. The control and driving 
electronics 57 may also include a suitable modulator, such as a phase shift keying (PSK) modulator, which conven- 
tionally produces driving signals for application to the transmitter 56. These driving signals can be used to apply ap- 
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propriate modulation to the mud siren of transmitter 56. It will be understood that alternative techniques can be employed 
for communicating logging information to the surface of the earth. 

The downhole subsystem 1 00 further includes acquisition and processor electronics 58. These electronics include 
a microprocessor (with associated memory clock circuitry, and interface circuitry), and processing circuitry. The ac- 

s quisition and processor electronics 58 are coupled to the measuring apparatus 200 and obtain measurement informa- 
tion therefrom. The acquisition and processor electronics 58 is capable of storing data from the measuring apparatus 
200, processing the data and storing the results, and coupling any desired portion of the information it contains to the 
transmitter control and driving electronics 57 for transmission to the surface by transmitter 56. A battery 53 may provide 
downhole power. As known in the art, a downhole generator (not shown) such as a so-called "mud turbine" powered 

10 by the drilling fluid, can also be utilized to provide power during drilling. If desired, the drilling equipment can optionally 
be a directional drilling system (not shown) which has, for example, a bottom hole assembly that includes a stabilizer, 
an offset (or "bent") sub, a mud motor that is driven by the flowing mud, and a near-bit stabilizer. The bent sub typically 
has an offset or bend angle of to 2 degrees. As is known in the art, when the bit is driven by the mud motor only 
(with the drill string stationary), the bit will deviate in a direction determined by the tool face direction in which the drill 

15 string and bent sub are oriented (so-called "sliding mode"). When it is desired to drill substantially straight, the drill 
string and the mud motor are both rotated at appropriate rates (so-called "rotating mode"). In this manner, directional 
drilling can be implemented with reasonable accuracy and without unduly frequent tripping of the drill string. 

Figure 2 is a block diagram of electronics in the measuring apparatus 200 and the acquisition and processor 
electronics 58. An oscillator 211 produces an electrical signal of frequency f (e.g. at 2 MHZ), which is amplified by 

20 amplifier 21 2 and applied, via electronically controlled switch 275, to one of the transmitting antennas T1 , T2, T3, T4 
or T5 through respective impedance matching circuits 221 . 222, 223, 224, and 225. The switch 275, under control of 
the processor 270, selects the transmitter T1 , T2, T3, T4 or T5 that is to be energized. The signal from receiver R1 is 
coupled, via impedance matching circuit 230 and preamplifier 231 . to one input of an electronically controlled switch 
250. Similarly, the signal from receiver R2 is coupled via impedance matching circuit 240 and preamplifier 241 to the 

25 other input of the switch 250. The switch 250, under control of processor 270, selects the (R1 ) or (R2) receiver output. 
The selected signal is amplified by amplifier 251 and subsequently translated to a lower frequency Af using a known 
heterodyne technique. A local oscillator 213 is in a phase locked loop (represented at 215) with the master oscillator 
211 . The local oscillator 213 has a frequency f+Af, where Af is typically a few kilohertz or less. A signal from the local 
oscillator 21 3 is mixed with the received signal by a mixer 255, and the mixer output is passed through a low pass filter 

30 257 which blocks signals at f and f+Af and passes the signal at Af . The signal at frequency Af contains the phase and 
amplitude information of the original signal at frequency f. The recognized advantage of this frequency translation is 
that it is easier to measure the phase and amplitude at kilohertz or lower frequencies than at megahertz frequencies. 
The signal at frequency Af is measured with a phase meter 261 and with an amplitude meter 262, and the results are 
input to the processor 270. The phase meter 261 may utilize a reference signal from the phase locked loop 215. The 

35 phase and amplitude meters may also include sample-and-hold circuits, to implement comparison of the signal from 
the respective receivers. Also, the processor 270 can compute the relative phase and amplitude (as well as phase 
average, as will be described) from the measurements it receives. Reference can be made to U.S. Patents 4, 1 85,238, 
4,899,112, and 5,594,343. 

As noted above, the processor 270 has associated memory, clocking, and interface circuitry (not shown), as is 

40 conventional. The processor 270 can implement storage of the measurements of phase and amplitude, processing of 
these measurements, storage of the results of the processing, and/or coupling of the measurements and/or processing 
results to the transmitter control and driving electronics 57 for transmission to the earth's surface. As is known in the 
art, a downhole clock can be utilized to keep track of time, which can subsequently be correlated with depth level by 
virtue of keeping a record of the drill string progress. (In general, depth level means longitudinal depth in the borehole.) 

45 The clock, which can typically be part of the system in conjunction with the processor 270, can be synchronized with 
the system before a trip downhole. Also, if desired, communication with the downhole subsystem, using mud pulse 
technology or other suitable communication means, can be utilized to convey timing synchronization and/or depth level 
information, as well as to communicate data in general. It will be understood that the invention can be utilized in con- 
junction with any suitable technique for keeping track of depth level. 

50 Consider the model space in the diagram of Fig. 3 in which the z direction is considered vertical and the x-y plane 

considered horizontal. Initially, assume that bedding planes are perpendicular to the z direction and let a v and a h be 
the conductivities in the vertical (z) direction and the horizontal direction (x-y plane), respectively. This is a so-called 
"transverse isotropic" (Tl) model. As described in Moran et al., "Effects Of Formation Anisotropy On Resistivity Logging 
Measurements", Geophysics , Vol. 44, No. 7, 1979, and summarized in the above-referenced U.S. Patent 5,329,448, 

55 the logging device coil fields can be considered as the superposition of magnetic dipoles 511 and 512 (Fig. 3) having 
respective magnetic moments. The magnetic dipole 511 is oriented in the vertical direction and the magnetic dipole 
512 is oriented in a direction in the horizontal plane (arbitrarily, the x direction, in this case). The following notation is 
introduced. 
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12 2 2 

r = a/jc + y + z 
f~2 2 72 2 

= vx + y + X z 
p = Vx + y 



(3) 



Then the Hertz vector and scalar potential for a magnetic dipole located atx = y = z = 0 and oriented in the x direction 
is given by: 
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(4) 



For a similar dipole oriented in the z direction: 



K x =0 



n = 0 
y 



71, = 



M z e h 



z 4n r 



_ A7 Z i^zf, 7 ^ 

<J) — J - e 

4k r 2 y ik h r J 
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(5) 



The electric and magnetic fields are given by: 



A 



1= = i<0[L Q c h ^ X 7C 



(6) 



7? = iCDjLl^^TC + ^<D 



(7) 



From these expressions for the electric and magnetic fields, the voltages at the receivers can be derived as being 
equal to m 2 nfi- H R • n where 7? is the direction along the axis of the tool, r R is the radius of the receiver coil, and R 
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is the number of turns on the receiver. For an induction or a propagation logging tool M- nR^-Ti where R T \s the radius 
of the transmitter coil, T is the number of turns on the transmitter coil, and I is the current. 

From equations (6) and (7) one can write the magnetic field at the receiver coil for a transmitter at the origin and 
the tool tilted at an angle G from z in the x-z plane. For simplicity, it can be assumed that the tool remains in the x-z 
plane. Then the components of the dipole moment of the transmitter are given by: 



M x = Msin 6 



(8) 



M z = M cos 9 



(9) 



Define H /y as the magnetic field in the i direction due to the component of the source in the j direction. This gives: 
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(13) 



Since the transmitter and receiver are aligned in the same direction the voltage in the receiver coil is given by: 



V= {m % r 2 R R)(cos Q(H ZX + H zz ) + sin Q( H 



(14) 



Once the voltage at each of the receiver coils due to each of the transmitter coils is determined, the total measurement 
can be determined by adding the voltages in the case of an induction tool ; or by taking the complex ratio of the voltages 
in the case of a propagation tool. For example, forthe propagation logging device of Figures 1 and 2, for each transmitter 
position, the absolute value of the voltage at each receiver can be obtained as the square root of the sum of squares 
of the real and imaginary parts of the complex voltage (equation (14)), and the ratio of the absolute values provides 
the attenuation, from which the attenuation-determined resistivity Rad is obtained. The phase for each receiver is 
obtained from the arc tangent of the ratio of the imaginary and real parts of the complex voltage. The phase shift is the 
difference in phase at the two receivers and the phase average is the average of the phases at the two receivers. The 
phase-shift-determined resistivity Rps and the phase-average-determined resistivity Rpa are respectively determined 
from the phase shift and the phase average. 

Shown in Figures 4 and 5 are examples of phase-shift (Rps), phase-average (Rpa), and attenuation (Rad) resistivity 
responses as a function of relative dip angle in a transverse isotropic (Tl) medium with relatively high (Rh=0.5 ohm-m 
and Rv=20 ohm-m) and relatively low (Rh=0.5 ohm-m and Rv=2 ohm-m) anisotropy ratio, respectively. These anisot- 
ropy response examples can be used to illustrate the following characteristics that can be expected for an anisotropy 
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inversion: 

1. Rps, Rpa, Rad hardly respond to Rv at low relative dip angle (Dip < 40 degree). At low relative dip angle, the 
values of these resistivities are very close to Rh regardless of whether the Rv value is 2 or 20 ohm-m. Therefore, 

s it can be expected that the inversion will have a poor resolution of Rv at low dip angle. 

2. At higher relative dip angle (dip>40 degree) ; the anisotropy response increases with dip angle. The Rps, Rpa, 
and Rad values may exceed Rv at high dip angle. At lower anisotropy ratio, the spreading of the Rps ; Rpa, and 
Rad as a function of transmitter-to-receivers spacing (TR) is orderly, namely, a larger TR measurement will read 
higher resistivity than a smaller TR measurement. However, at relatively high anisotropy ratio, the spreading of 

10 the Rps, Rpa, and Rad as a function of TR is not orderly. This is illustrated for example in Figure 4 where the Rps 

34 inches may be lower then Rps 28 inches, etc. This non-orderly structure precludes use of a simple interpolation 
type of inversion algorithm. A nonlinear least square routine, such as a Levenberg-Marquardt algorithm, together 
with an extensive initial search routine can be effectively used to handle this difficulty. 

3. For cases of high anisotropy ratio at high relative dip angle such as that shown in Figure 4, most of the Rps and 
15 Rad values may go out of the resistivity transform range due to negative phase shift and below geometric spreading 

attenuation responses. In these cases, it is necessary to include Rpa in the inversion algorithm to enhance the 
robustness of the algorithm. For example, in the case of Figure 4 above 85 degrees there are only two valid 
measurements in Rps and Rad (Rps 34 inches and Rad 10 inches) which are insufficient to invert for three pa- 
rameters, Rh, Rv. and Dip without using Rpa measurements. 

20 

Referring to Figure 6, there is shown a flow diagram of a routine which, in conjunction with the routines referred 
to therein, can be used in programming the processors, such as processors 85 and 270 of Figures 1 and 2, respectively 
to implement operation in accordance with an embodiment of the invention. In the illustrative embodiment hereof, 
certain operations are performed and/or controlled by downhole or by uphole electronics, as the case may be, but it 

25 will be understood that the invention can be implemented with some of the operations performed at different alternative 
locations. The block610 represents the routine, described in further detail in conjunction with Figure 7, for implementing 
the borehole logging and storage of measurements. This includes measurement of amplitude and phase at each re- 
ceiver for each of the transmitter-to-receivers spacings. (As described in U.S. Patent 5,594,343, in addition to the five 
actual transmitter-to-receivers spacings, two more virtual spacings can be computed, so a total of seven spacings are 

30 available.) In the present embodiment the measured phase shift PS m , the measured phase average PA m , and the 
measured attenuation AD m , are computed for each of the transmitter-to-receivers spacings. 

The block 620 represents implementing of a setup routine wherein operator selected parameters can be input, if 
desired. 

For example, weighting values, treated below, can be selected. Also, if desired, screening can be performed, either 

35 manually or automatically, to predetermine (e.g. with a preprocessing routine) conditions under which the technique 
hereof would not be used or would require modification. An example of such a condition would be when there is no 
meaningful anisotropy, as may be indicated by substantial correspondence in measurements obtained from all of the 
various transmitter-to-receivers spacings of the logging apparatus. 

The block 630 represents the computation, for each transmitter-to-receivers spacing i, and over the depth range 

40 of interest, of the resistivity determined from the measured phase shift, Rps_m(i), resistivity determined from the meas- 
ured phase average, Rpa_m(i), and resistivity determined from measured attenuation Rad_m(i). The block 635 rep- 
resents initializing to a first depth level in the depth range of interest. The block 635 is then entered, this block repre- 
senting a routine, described in conjunction with the flow diagram of Figure 8, for implementing an initial search to 
determine candidates for initial model values of Rh, Rv, and Dip. This search is global in nature and is used to provide 

45 candidate model values that are close enough to the correct final solution so that the next more refined process of 
inversion is likely to converge to the correct final solution rather than to another solution made possible by the complex 
(and not always orderly) nature of the relationships between the measured resistivities at various spacings and the 
anisotropic formation dip angle. The block 650 is then entered, this block representing a routine, described in conjunc- 
tion with Figure 9, for implementing the inversion to determine model values of Rh, Rv and Dip. The determined values 

50 are read out and recorded as represented by the block 660. The decision block 670 is then entered, this block repre- 
senting inquiry as to whether the last depth level to be processed in the present pass has been reached. If not, the 
depth level index is incremented (block 675), block 640 is reentered, and the loop 678 continues until all desired depth 
levels of the depth range have been processed. Readout can be, for example, to another storage or recording medium, 
or to a display (not shown). (Although all the processing is shown as being performed at the borehole site in the 

55 illustrated embodiment, it will be understood that, if appropriate, the processing, or a portion thereof, can be performed 
remote from the borehole site, such as by a communications link.) 

Referring to Figure 7, there is shown a flow diagram of a routine, represented generally by the block 610 of Figure 
6, for controlling a processor (e.g. 270 of Figure 2) of the embodiment of Figures 1 and 2, to implement borehole logging 
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and store measurements. (The same type of routine can be used for any number of transmitters.) The blocks 710, 
720, 730 ; and 740, and the loop 745, operate to sequentially cause energizing of each transmitter (T1 through T5, in 
this case), and the controlling of the measuring, at each receiver, of the relative phase and the amplitude of the received 
signals associated with each energized transmitter. Reference can be made to the above referenced U.S. Patent 

s 5,594,343 for further details. It will be understood also that other techniques can be employed for implementing the 
measurements. The block 710 represents initializing to the first transmitter and the block 720 represents activating 
the transmitter (by control of switch 275 of Figure 2) and storing the relative phase and the amplitude measured at 
each of the receivers of the receiver pair R1 and R2 (see Figure 2). The decision block 730 represents the determination 
of whether the last transmitter has been activated. If not, the block 740 is entered, the transmitter index is incremented, 

10 and the block 720 is reentered. The loop 745 continues until all measurements have been made and stored for the 
present depth level. 

The blocks 770, 780 and 790 are used to respectively compute, for each transmitter-to-receivers spacing, the 
attenuation, AD, the phase shift, PS, and the phase average, PA, from the relative phase and amplitude measurements 
at near receiver R1 (respectively designated 0-, and A-,) and the relative phase and amplitude measurements at far 
15 receiver R2 (respectively designated cf> 2 and A 2 ). In the present embodiment, the quantities AD, PS, and PA are com- 
puted from the following relationships: 

AD = 20log 10 [A.,/A 2 ] 

20 

PS = 4> 2 -0 1 



25 PA=W 2 +^y2 

The block 760 represents initializing to the first transmitter-to-receivers spacing, whereupon AD, PS and PA are com- 
puted, as indicated. The decision block 795, increment block 798, and loop 797 are then continued until all spacings 
have been processed. Alternatively the blocks 770, 780, and 790 could be under the block 720 in the loop 745, for 

30 computation of AD, PS, and PA directly after each set of measurements for a particular spacing is taken. 

In an embodiment hereof, the inversion routine utilized is a Levenberg-Marquardt type of algorithm. Starting from 
initial guess values (determined in a manner to be described), a series of iterations is used in following a path of 
steepest descent to a local minimum. The inversion technique is described in conjunction with the routine of Figure 9, 
but is first summarized as follows: 

35 Rps_m(i), Rpa_m(i), Rad_m(i), are the measured phase-shift, phase-average, and attenuation resistivities at the 

i-th mixed borehole compensated spacing TR(i), i=1 ,2,...,ntr. In an embodiment hereof, TR(i) = [10, 13, 16, 22, 28, 31 , 
34] inches, ntr=7. 

Define Rps(Rh, Rv, Dip, TR(i)), Rpa(Rh, Rv, Dip, TR(i)), and Rad(Rh, Rv, Dip, TR(i)) as the model computed phase- 
shift, phase-average, and attenuation resistivities at the i-th mixed borehole compensation spacing TR(i) : i=1,2 ; ..., ntr. 
40 Since the basic measurements of phase-shift, phase-average, and attenuation are all responding to conductivity 

it is advantageous to match the measurements with model values in conductivity space such that the over-all error in 
matching will not be dominated by a few high resistivity data points. Therefore, define a set of conductivity terms 
corresponding to the above resistivity terms: 

45 

Cps_m(i) = 1/Rps_m(i) 



Cpa_m(i) = 1/Rpa_m(i) 

50 

Cad_m(i) 1/Rad_m(i) 



Cps(Rh, Rv, Dip, Tr(i)) = 1/Rps(Rh, Rv, Dip, TR(i)) 
Cpa(Rh, Rv, Dip, Tr(i)) = 1/Rpa(Rh, Rv, Dip, TR(i)) 
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Cad(Rh, Rv, Dip, Tr(i)) = 1/Rad(Rh, Rv 3 Dip, TR(i)). 



The error spanned between the measured data and modeled data is given by: 



ntr 

Err <Rh, Rv, Dip) = S {Wps(i) * [Cps(Rh, Rv, Dip, TR ( i ) ) - 

i=l 

Cps_m ( i ) ] 2 + 
Wpa(i) * [Cpa(Rh, Rv, Dip, TR < i ) ) - 

Cpa_m ( i ) ] 2 + 
Wad(i) * [Cad(Rh, Rv, Dip, TR ( i ) ) - 

Cad_m(i)] 2 } 

(15) 



Here, Wps(i), Wpa(i), Wad(i), i=1,...,ntr are optional weighting coefficients for phase-shift, phase-average, and atten- 
uation components, respectively. If, for example, the weighting coefficients of any one component are set to zero, the 
inversion results can then be made independent of that component's data. In addition to assigning a relative importance 
to the data components, the weighting coefficients can also be used as a switch to shut off any combination of the 
components. For example, setting Wps(i)=1, Wpa(i)=1, Wad(i) =0, i=1,..., ntr, would mean only using the phase-shift 
and phase-average to do the inversion. 

The routine hereof (e.g. a Levenberg-Marquardt algorithm in the present embodiment) will determine the values 
of inverted model parameters Rh_inv, Rv_inv and Dip_inv such that Err(Rh_inv, Rv_inv, Dip_inv) achieves a local 
minimum. The routine requires an initial guess of the model parameters. 

(As discussed elsewhere herein the general accuracy of the initial guess or guesses is quite important.) Starting 
from the initial guess values, through a series of iteration steps, the algorithm will follow the path of steepest descent 
to a local minimum. During the iteration steps, there may be times when the algorithm may determine that the next 
model parameters to try are negative Rh and/or Rv values. It is therefore necessary to impose a stiff exponentially 
increased Err value when the guessed Rh and/or Rv values fall below certain values to prevent this condition from 
occurring. An example of the modified error function, Err_m, is given as: If Rh > 0.05 and Rv > 0.05 then 



ntr 

Err_m (Rh, Rv, Dip) = 2 {Wps(i) * [Cps(Rh, Rv, Dip, TR(i) ) - 

i=l 

Cps_m ( i ) ] 2 + 
Wpa(i) * [Cpa(Rh, Rv, Dip, TR ( i ) ) - 

Cpa_m ( i ) ] 2 + 
Wad(i) * [Cad(Rh, Rv, Dip, TR { i ) ) - 
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Cad_m ( i ) ] 2 } 

s otherwise 



Let R = min(Rh, Rv) 
if R > 0,02 then 

fexp = 0.5/R 

otherwise 
fexp =25.0 
and 

ntr 

Err_m (Rh, Rv, Dip) = S {Wps(i) * [exp(fexp)] 2 + 

i=l 

Wpa(i) * [exp(fexp)] 2 + 
Wad(i) * [exp(fexp) ] 2 } 



Figure 8 is a flow diagram of the routine represented by the block 640 of Figure 6A for implementing the initial 
search to obtain candidates as initial model values to be used in the subsequent inversion routine. A coarse grid of 
30 values is set up for the model parameters Rh, Rv and Dip, as represented by the block 810. In the present embodiment, 
the grid points for Rh and Rv are denoted by R_grid which is a 22-point vector given as: 

R_grid = [0.1, 0.15, 0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 3.0, 5.0, 7.0, 10.0, 15, 20, 30, 50, 70 ; 100, 150, 200, 300] ohm-m, 

and the grid points for the relative dip angle are denoted by Dip_grid which is a 36-point vector given as: Dip_grid = 

[0, 10, 20, 30, 40, 50, 52 ; 54, 56, 58, 60, 62, 64, 66, 68, 70, 71 , 72, 73, 74, 75, 76, 77, 78, 79, 80, 81 , 82, 83, 84, 85, 
35 86, 87, 88, 89, 90] degrees. The previously obtained phase-shift, phase average, and attenuation resistivities, Rps_m 

(i), Rpa_m(i) and Rad_m(i), i = 1 ntr are retrieved, as represented by the block 820. Optional weighting factors for 

the resistivities, Wps(i), Wpa(i), and Wad(i), can also be retrieved. 

Based on the Rps_m(i) values, the present routine determines the range of Rh and Rv and Dip to be searched 

over the coarse grid. Then the routine finds three candidates (Rh_ig(m), Rv_ig(m), Dip_ig(m), m=1, 2, 3) over the 
40 search range which have the smallest values of Err_m(Rh_ig(m), Rv_ig(m), Dip_ig(m)) as defined in equation (16). 

The rules, for the present embodiment, for determining the search range for the model parameters are as follows: 

Let Rps_m_max, and Rps_m_min be the maximum and minimum values of Rps_m(i), i=1,...,ntr, and let 
Rps_m_ratio = Rps_m_max/Rps_m_min. 
45 The boundary points for the search range for the model parameters Rh, Rv, Dip, denoted by Rh_start, Rh_end, 

Rv_start, Rv_end, Dip_start, Dip_end, are defined by the following: 
Rh_end = Rps_m_min 

if Rps_m_ratio > 2 or Rps_m_min = 1000 and Rps_m_max = 1000, then 
Rh_start = max (Rps_m_m in/500, 0.1) otherwise 
50 Rh_start = max (Rps_m_min/10, 0.1), and; if Rps_m_min = 1000 and Rps_m_max = 1000, then 

Rv_start = 20 
Rv_end = 1000 otherwise 
Rv_start = Rps_m_min 

Rv_end = min (Rps_m_min*10, 1000), and; Dip starts 30 
55 Dip_end = 90. 

Referring again to Figure 8, the block 825 represents determination of the maximum and minimum values of Rps_m 
(i) and the block 830 represents the determination of Rps_m_ratio of the maximum to the minimum. The end of the 



25 



(16) 
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range for Rh, designated Rh_end, is set equal to Rps_m_min (block 835). Determination is then made (decision block 
837) as to whether Rps_m_min and Rps_m_max are both equal to 1000 ohm-m. (In the present embodiment, 1000 
ohm-m is the highest resistivity used, so when Rps_m_min is 1000 ohm-m, Rps_m_max will also be 1000 ohm-m.) If 
so, Rh_start (the start of the Rh range) is set to the greater of Rps_m_min/500 or 0.1 ohm-m (block 839), Rv_start is 
set to 20 ohm-m (block 841 ), and Rv_end is set to 1 000 ohm_m (block 843). If not, the Rps_m_ratio is checked (decision 
block 838) to determine if it is greater than 2. If so, the block 839 is entered and the values or Rh_start, Rv_start, and 
Rv_end will be determined by the blocks 839, 841 and 843, as above. Otherwise, Rh_start is set to the greater of 
Rps_m_min/10 or 0.1 ohm-m (block 845), Rv_start is set equal to Rps_m_min (block 847), and Rv_end is set to the 
lesser of Rps_m_min*10 or 1000 ohm-m (block 849). Also, Dip_start is set to 30 degrees and Dip_end is set to 90 
degrees (block 851). Next, Rh(i), Rv(j) and Dip(k) are initialized within the selected ranges (block 860) and EL2(i,j ; k) 
is computed for all i,j.k in the range, in accordance with equation (16), (block 862, in conjunction with nested loops 865 
(for i, including blocks 863 and 864), 870 (for j, including blocks 867, 868 and 869) and 877 (for k, including blocks 
872, 874 and 875)). The three smallest EL2(i,j, k) are then determined (block 880). Next, the candidates for initial guess 
are determined as Rh(i), Rv(j), and Dip(k) that produced the three smallest EL(i,j ; k) (block 890). 

In most of the anisotropy conditions created in the natural deposition process, Rh is smaller than Rv. Therefore, 
the initial guess algorithm described here is tailored to find such a solution, although theoretical solutions with Rh > 
Rv frequently exist. 

Let Rh(i), Rv(j), Dip(k) be the model parameter grid point values which fall within the boundary points of the model 
parameters defined above, then 

Rh_start < Rh(i) < Rh-end, i-1 ,2,..., nrh 



Rv_start < Rv(i) < Rv-end, j-1 ,2,...,nrv 



Dip_start < Dip(k) < Dip_end, k=1 ,2,...ndip 

Let EL2(i,j,k) be the Euclidian L2 norm between the measured and the model data at all the grid points within the 
search range, then we have 



ntr 

EL2 (i, j ,k) = S {Wps(n) * [Cps(Rh(i), Rv(j), Dip(k), TR (n) ) - 

n=l 

Cps_m(n} ] 2 + 

Wpa(n) * [Cpa(Rh(i), Rv(j), Dip(k), TR(n)) - 

Cpa_m ( n ) ] 2 + 

Wad(n) * [Cad(Rh(i), Rv(j), Dip(k>, TR (n) ) - 

Cad_m(n) ] 2 } 

i=l , 2 , . . . , nrh, j=l, 2 , . . . , nrv, k=l , 2 , . . . , ndip 

(17) 



The candidates for initial guess of the model parameters are those that yield the three smallest EL2(i,j,k): 
EL2 (Rh_ig(1), Rv_ig(1), Dip_ig(1))< EL2(Rh_ig(2) ; Rv_ig(2), Dip_ig(2)) < 
EL2 (Rhjg(3), Rv_ig(3), Dip_ig(3))< EL2(i,j : k). 

for all i,j,k such that [Rh(i), Rv(j), Dip(k)] is not equal to [Rh_ig(m), Rv_ig(m), Dip_ig(m)], m=1 ,2,3. 

Figure 9 is a flow diagram of the routine, represented by the block 650 of Figure 6B, for implementing inversion 
to determine the inverted model values of Rh, Rv, and Dip. The initial guess index, m ; is initialized (block 905) and the 
initial guess model parameters are started at Rh_ig(m), Rv_ig(m), and Dip_ig(m) (block 910). Next, as represented by 
the block 920, the error value Err (Rh, Rv, Dip) is computed in accordance with equation (16). As previously noted, the 
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algorithm used for this part of the processing in the present embodiment, is a Levenberg-Marquardt algorithm, which 
is a well known inversion algorithm. Reference can be made, for example, to Numerical Analysis - The Levenberg- 
Marquardt Algorithm, Implementation and Theory , by Jorge J. More, G.A. Watson, Editor, Springer- Verlag, 1977, and 
to U.S. Patent 5,329,448 which utilizes a similar general technique of inversion by following a path of steepest descent 

s to a minimum. The decision block 930 represents determination of whether a predetermined error criterion has been 
met or, if desired, whether a maximum number of iterations have been performed. If not, the block 932 is entered, this 
block representing the computation of the direction and step size of the increment. The model parameters are then 
incremented (block 934), the block 920 is re-entered, and the loop 935 continues until the inquiry of block 930 is 
answered in the affirmative. The inverted model parameters for the current m are determined as Rh_lm(m), Rv_lm(m) 

10 and Dip_lm(m), along with the L(2)norm as a measure of the fit (block 940). Inquiry is then made (decision block 950) 
as to whether the maximum for index m (3, for this embodiment) has been reached. If not m is incremented (block 
955), block 910 is re-entered, and the loop 957 is continued until the inversion has been implemented starting from 
each of the three initial guess model value parameter sets. When this has been done, the smallest L2_lm(m) is selected 
(block 960) and the inverted solution model value parameters Rh_inv, Rv_inv, and Dip_inv are selected as the ones 

15 with m corresponding with the m of the smallest L2_lm(m) (block 970). 

In the above processing there can be implicit assumptions regarding both the horizontal and vertical dielectric 
constants, which can be expressed as respective functions of horizontal and vertical resistivity, for example 

20 eh = 108.5 • Rh " 0 35 (18) 

ev= 108.5 *Rv (19) 

25 

These assumptions provide additional equations such that one can invert for Rh, Rv, and relative dip angle (Dip) from 
a minimum of 3 measurements out of the available Rps, Rad, and Rpa measurements at various spacings. 

Equations (18) and (1 9) are derived empirically based on test data from a limited number of rock samples. Although 
these assumptions are reasonably good, if the actual values of eh and ev of a given formation deviate significantly 
from those predicted by equations (18) and (19), substantial error in Rv and Dip may result. 

30 

In an embodiment hereof, inversion can be implemented for Rh and eh without invoking any dielectric assumption. 
Rv and ev can then be calculated if the Dip angle is known. If the Dip angle is unknown, but is within the range of 80° 
< Dip < 100° (a condition commonly met for most horizontal wells) and the degree of anisotropy is not very large, Rv 
and ev and Dip can be calculated if one makes a dielectric assumption, such as equation (1 9), on the vertical component. 
This alternative method allows the flexibility of delaying the making of any assumption until after the inversion. 

35 

For a general case, the signals at the receiver coils (induced voltage, phase shift, and attenuation) can be expressed 
as functions of the five unknown parameters, ah, eh, ov, ev, and G in the following functional form (see e.g. "Processing 
and Modeling 2-Mhz Resistivity Tools in Dipping, Laminated, Anisotropic Formations" by Martin G. Luling, Richard A. 
Rosthal, and Frank Shray, SPWLA35th Annual Logging Symposium, June 19-22, 1994, Paper QQ; and "A New Method 
to Determine Horizontal-Resistivity in Anisotropic Formations Without Prior Knowledge of Relative Dip" by T Hagiwara, 
SPWLA37th Annual Logging Symposium, June 16-19, 1996, Paper Q] : 

V = V(ah, eh, ov, ev, 6, TR) = V (Kh, p, TR) (20) 

45 

Kh = (icou.Cn) 172 (21) 

P = [cos 2 G + (Cv/Ch)sin 2 6] 1/2 (22) 

50 

Ch= ah - icoeh (23) 

55 Cv=av-icoev (24) 

The symbols used in equations (20)-(24) are defined as the follows: 
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V -the induced the voltage at the receiver coil 

TR -the transmitter-to-receivers distance 

Kh -the horizontal wave number 

p -the anisotropy factor 

Ch -the complex horizontal conductivity 

Cv -the complex vertical conductivity 

i -the imaginary operator 

co -the angular frequency (=27if, f is frequency) 

jlx -the magnetic permeability 

With four or more measurements, one can invert for the two complex numbers, Kh and p. From the real and 
imaginary part of Kh, one can solve for Ch and then oh and eh using equations (21) and (23). From the real and 
imaginary part of p and Ch, there are not a sufficient number of equations to solve for the remaining three unknown 
parameters, ay ev, and 6. However, without implementing any further inversion, the following options can be used to 
obtain further results: 

(a) Use D & I (direction and inclination) data, such as from a measuring-while-drilling tool and local knowledge of 
the structural dip of the formation to estimate the relative dip angle 9 and treat it as a known input. Then, using 
equations (22) and (24), one can solve for av and ev from the given values of p and 0. 

(b) Without solving for 6, one can still estimate fairly accurately the values of gv and ev if the degree of anisotropy 
is not excessive and the value of 6 is within the range of 80° < 6 <100°. This range of 9 occurs commonly for the 
condition of a horizontal well through horizontal beds. From equation (22), one can obtain, 

A + iB = (p 2 - 1 )Ch = (Cv-Ch)sin 2 G (25) 
A= (ov-oh)sin 2 9 (26) 

B= -co(ev-eh)sin 2 9 (27) 

From the values of p and Ch, one can solve for the values of A and B which, in turn allows one to solve for av 
and ev, respectively, through equations (26) and (27). In these two equations, the dependence on 6 is through the 
term sin 2 9, which changes very slowly for 9 near 90°, i.e. for 9 in the range of 80° < 9 < 100°. Within this range, 
one can obtain a good approximation of the values of av and ev by assigning 9 = 90°: 

av = ah + A (28) 



ev = eh - B/co (29) 

These approximated values av and ev are always higher than the true values of av and ev, respectively. It can be 
shown that the fractional errors of the approximated vertical conductivity and dielectric constant can be written as 
a function of 9 and their respective ratios of horizontal to vertical values: 

(av - av)/av = (1 - ah/av)(sin 2 9 - 1) (30) 

(ev - ev)/ev= (1 - eh/ev) (sin 2 9 - 1) (31) 

The fractional error for the approximated vertical conductivity and dielectric constant are plotted on Figure 10. It 
can be seen that the fractional errors for the approximated vertical conductivity and dielectric constant are indeed 
small for 80° < 9 < 100° and moderate horizontal to vertical ratio (<10). 

(c) Assuming a relationship between av and ev, such as equation (1 9). Such a dielectric relationship, together with 
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equation (22), provides sufficient equations to solve for the remaining three unknown parameters, ov, ev, and 0. It 
is more advantageous to invoke a dielectric assumption after the inversion for Kh and p is finished instead of 
committing to specific assumptions first then inverting for the remaining parameters. After obtaining Kh and p, one 
can repeatedly try different assumptions, if necessary, to determine the most consistent one without incurring the 
cost of repeat inversions. 

The approach of the inversion of the present embodiment is similar to that of the earlier embodiment, but in this 
case with a formation model which contains four unknown parameters, Khr, Khi, pr, and pi. Here, Khr and Khi are the 
real and imaginary parts of the anisotropy factor p. In this case PS(Khr, Khi, pr, pi, TR(i)), PA(Khr, Khi, pr, pi TR(i)) and 
AD(Khr, Khi, Pr, pi, TR(i)) are the model computed phase-shift, phase -average, and attenuation at the i-th mixed bore- 
hole compensation spacing TR(i), i=1,2,...,ntr. In this more general case, the error spanned between the measured 
data and modeled data is given by: 



ntr 

Err (Khr, Khi, pr, pi) = 2 { Wps(i) * [ PS { Khr , Khi, Pr, pi, 

i=l TR(i) ) - PS_m(i) ] 



2 



Wpa(i) * [PA(Khr, Khi, pr, pi, 

TR ( i ) ) - PA_m ( i ) ] 2 + 

Wad(i) * [AD (Khr , Khi, Pr, Pi, 

TR(i) ) - AD_m(i) ] 2 } (32) 



For the above exemplary logging tool there are a total of twenty one measurements (PS_m(i), PA_m(i) and AD_m 
(i)) available for the inversion process. Since there are only four unknown parameters in the model, a minimum of four 
measurements are needed for the inversion. If all the measurements are of the same high quality (low noise and fit 
the Tl model), due to the benefit of additional averaging, the more measurements used in the inversion the better the 
results. However, if some measurements are corrupted more than others due to noise or deviation from the Tl model, 
then inclusion of the more corrupted data in the inversion may lead to a lower quality result. Therefore, one should 
preferably choose only the highest quality measurements (at least four of them) to be included in the inversion process. 

In this case, the initial search algorithm searches for the initial guess values for model parameters Khr, Khi, pr, pi 
through searching over a coarse grid values of Rh, Rv, eh, ev and 9. The grid point values for eh and ev are determined 
by the grid point values of Rh and Rv, respectively, through equations (18) and (19). A 10-point grid is set up for the 
relative dip angle 9. The values of the Rh, Rv, and 6 grid for this case can be: 

Rh(i)=10 (-1 0 + i-0.04) ohm-m, i=1,2.,„100 
Rv(j)=10 (-i.o + i-0.04) ohm-m, j=1,2.,,,100 

Dip(k)= [0, 50, 55, 60, 65, 70, 75, 80, 85, 90] degrees, for k=1,2,...,10 
Also for this case the model parameters Khr, Khi, pr, pi at the grid points can be written as: 



Khr_g(i) = Khr(Rh_g(i)) 
Khi_g(i) = Khi(eh(Rh_g(i)) 

pr_g(i,j,k) = pr(Rh_g(i), Rv_g(j), eh(Rh_g(i)), ev(Rv_g(j), 
Dip_g(k)) 

pi_g(i,j,k) = Pi(Rh_g(i), Rv_g(j), eh(Rh_g(i», ev(Rv_g(j), 
Dip_g(k)) 

i=1 ,2,...nrh, j=1 ,2,...,nrv, k=1 ,2,...,ndip 



and the Euclidian L2 norm between the measured and the model data at all the grid points within the search range is: 
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ntr 

EL2(i,j,k) = £ { Wps(n) * [ PS (Khr_g ( i ) , Khi_g(i), 
pr_g(i,j ,k) , n=l 

Pi_g(i, j ,k) , TR (n) ) -PS_m (n) ] 2 
+ Wpa(n)*[PA(Khr_g(i) , Khi_g(i), pr_g(i,j,k), 

pi_g(i,j,k) , TR ( n ) ) - PA_m ( n ) ] 2 

+ Wad(n) * [AD(Khr_g(i) , Khi_g(i), Pr_g (i , j , k) , 

Pi_g(i,j,k) , TR (n) ) - 

AD_m(n) ] 2 } 

i = l , 2 , . . . nrh, j=l, 2 , . . . , nrv, k=l , 2 , . . . , ndip 



Again, the candidates for initial guess of the model parameters are those that yield the three smallest EL2(i,j,k). 
The candidates for the initial guess of the four model parameters are given as: 



Khr_g(i1), Khi_g(i1), Pr_g(i1,j1,k1), PLg(i1J1,k1) 



Khr_g(i2), Khi_g(i2), |3r_g(i2 ; j2 ) k2) ) pi_g(i2J2,k2) 



Khr_g(i3), Khi_g(i3), pr_g(i3J3,k3), pi_g(i3J3,k3) 

The final inverted solution is chosen as the one that has the smallest L2_lm value: 

Khr_inv = Khrjm(m), 
Khi_inv = Khi_lm(m), 
pr_inv = prjm(m), 
pi_inv = pi_lm(m), 
L2_inv = L2Jm(m), 

where m is the index such that L2_lm(m) = min(L2_lm(n), n=1 ; 2,3). 
L2_inv can be used as a quality indicator for the inversion. 

After obtaining the values of Khr, Khi, pr, and pi, Rh and eh can be obtained directly from Khrand Khi using equation 
(21) and (23). Rv, ev, and, for some cases, 6, can be obtained according to the three options outlined above. 

(a) Assuming 0 is known: 

Rv and ev can be computed using equations (21) and (23). 

(b) When 9 is unknown but is within the range 80° < 0 < 100° and moderate horizontal to vertical ratio (<10): Rv 
and ev can be approximated using equations (23), (25), (26), (27), (28), and (29). 

(c) When a dielectric assumption on the vertical component, such as equation (1 9) can be made: 
Rv, ev, and 0 can be computed using equations (1 9) and (22). 
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Claims 



1. A method for determining the horizontal resistivity, vertical resistivity, and dip of formations surrounding a borehole, 
comprising the steps of: 

5 

(a) suspending a logging device in the borehole; 

(b) transmitting electromagnetic energy from a transmitter location on the logging device, and receiving the 
transmitted electromagnetic energy at receiver locations on the logging device for a first transmitter-to- receiv- 
ers spacing associated with said transmitter and receiver locations; 

10 (c) determining, from the received electromagnetic energy, measurement characteristics associated with said 

first transmitter-to-receivers spacing; 

(d) repeating steps (b) and (c) for a plurality of further transmitter-to-receivers spacings to obtain measurement 
characteristics for said plurality of further transmitter-to-receivers spacings; 

(e) generating a model of anisotropic formations that has horizontal resistivity Rh, vertical resistivity Rv, and 
15 a Dip angle with respect to a borehole reference; 

(f) selecting initial model values of Rh, Rv and Dip; 

(g) computing an error value from the differences, for each of a plurality of said transmitter-to-receivers spac- 
ings, between said measurement characteristics and model measurement characteristics obtained from said 
model using the model values of Rh, Rv, and Dip; 

20 (h) modifying said model values of Rh, Rv and Dip; 

(i) repeating steps (g) and (h) to minimize said error value: and 

(j) outputting the ultimately modified model values of Rh, Rv, and Dip. 

2. The method as defined by claim 1 , wherein said measurement characteristics comprise phase-shift, phase-aver- 
25 age, and attenuation. 

3. The method as defined by claim 1, wherein said measurement characteristics comprise resistivity derived from 
phase-shift, resistivity derived from phase-average, and resistivity derived from attenuation. 

30 4. The method as defined by claim 1 , wherein said model of anisotropic formations is a transverse isotropic model. 

5. The method as defined by claim 2, further comprising repeating said method at a multiplicity of depth level positions 
of said logging device in said borehole to obtain a log of at least one of Rh, Rv and Dip. 

35 6. The method as defined by claim 2, further comprising repeating said method at a multiplicity of depth level positions 
of said logging device in said borehole to obtain logs of Rh, Rv and Dip. 

7. The method as defined by claim 1, wherein said step of suspending a logging device in the borehole comprises 
coupling a logging-while-drilling device in a drill string in the borehole. 

40 

8. The method as defined by claim 1 , wherein the step of selecting initial model values of Rh, Rv, and Dip comprises: 

producing discretized model values of Rh, Rv and Dip; 

comparing said measurement characteristics with model measurement characteristics obtained using a mul- 
45 tiplicity of combinations of discretized model values; and 

selecting, as said initial model values, discretized model values based on the comparisons of said comparing 
step. 



9. The method as defined by claim 2, wherein the step of selecting initial model values of Rh, Rv, and Dip comprises: 

50 

producing discretized model values of Rh, Rv and Dip; 

comparing said measurement characteristics with model measurement characteristics obtained using a mul- 
tiplicity of combinations of discretized model values; and 

selecting, as said initial model values, discretized model values based on the comparisons of said comparing 
55 step. 



10. A method for determining properties of anisotropic formations surrounding a borehole, comprising the steps of: 
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(a) suspending a logging device in the borehole; 

(b) transmitting electromagnetic energy from a transmitter location on the logging device, and receiving the 
transmitted electromagnetic energy at receiver locations on the logging device for a first transmitter-to-receiv- 
ers spacing associated with said transmitter and receiver locations; 

s (c) determining, from the received electromagnetic energy, measurement characteristics associated with said 

first transmitter-to-receivers spacing; 

(d) repeating steps (b) and (c) for a plurality of further transmitter-to-receivers spacings to obtain measurement 
characteristics for said plurality of further transmitter-to-receivers spacings; 

(e) generating a model of anisotropic formations that has model anisotropy properties; 
10 (f) selecting initial model values of said anisotropy properties; 

(g) computing an error value from the differences, for each of a plurality of said transmitter-to-receivers spac- 
ings, between said measurement characteristics and model measurement characteristics obtained from said 
model using the model values of said anisotropy properties; 

(h) modifying said model values of said anisotropy properties; 
15 (i) repeating steps (g) and (h) to minimize said error value; and 

(j) outputting the ultimately modified model values of said anisotropy properties: 

said step of selecting initial model values of said anisotropy properties including the steps of: 
(k) producing discretized model values of said anisotropy properties; 

(I) comparing said measurement characteristics with model measurement characteristics obtained using a 
multiplicity of combinations of discretized model values; and 

(m) selecting, as said initial model values, discretized model values based on the comparisons of said com- 
paring step. 

11. The method as defined by claim 10 ; wherein said measurement characteristics comprise phase-shift, phase-av- 
erage, and attenuation. 

12. The method as defined by claim 10, wherein said measurement characteristics comprise resistivity derived from 
30 phase-shift, resistivity derived from phase-average, and resistivity derived from attenuation. 

13. The method as defined by claim 10, wherein said model of anisotropic formations is a transverse isotropic model. 

14. The method as defined by claim 10, further comprising repeating said method at a multiplicity of depth level posi- 
35 tions of said logging device in said borehole to obtain a log of at least one of horizontal resistivity Rh 3 vertical 

resistivity Rv and Dip angle with respect to a borehole reference. 

15. The method as defined by claim 10, wherein said step of suspending a logging device in the borehole comprises 
coupling a logging-while-drilling device in a drill string in the borehole. 

40 

16. The method as defined by claim 10, wherein said anisotropy properties are Khr, Khi, pr and pi, where Khr and Khi 
are the real and imaginary parts of the horizontal wave number Kh, and pr and pi are the real and imaginary parts 
of the anisotropy factor p, and where Kh is a function of horizontal conductivity and horizontal dielectric constant 
and p is a function of vertical conductivity and vertical dielectric constant. 

45 

17. The method as defined by claim 11 , wherein said anisotropy properties are Khr, Khi, pr and pi, where Khr and Khi 
are the real and imaginary parts of the horizontal wave number Kh, and pr and pi are the real and imaginary parts 
of the anisotropy factor p, and where Kh is a function of horizontal conductivity and horizontal dielectric constant 
and p is a function of vertical conductivity and vertical dielectric constant. 

50 

18. The method as defined by claim 1 3, wherein said anisotropy properties are Khr, Khi, pr and pi, where Khr and Khi 
are the real and imaginary parts of the horizontal wave number Kh, and pr and pi are the real and imaginary parts 
of the anisotropy factor p, and where Kh is a function of horizontal conductivity and horizontal dielectric constant 
and p is a function of vertical conductivity and vertical dielectric constant. 

55 

19. Apparatus for determining the horizontal resistivity, vertical resistivity, and dip of formations surrounding a borehole, 
comprising: 
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(a) a logging device suspendible in the borehole; 

(b) means for transmitting electromagnetic energy from a transmitter location on the logging device, and re- 
ceiving the transmitted electromagnetic energy at receiver locations on the logging device for a first transmitter- 
to-receivers spacing associated with said transmitter and receiver locations; 

(c) means for determining, from the received electromagnetic energy measurement characteristics associated 
with said first transmitter-to-receivers spacing; 

(d) means for repeating operation of elements (b) and (c) for a plurality of further transmitter-to-receivers 
spacings to obtain measurement characteristics for said plurality of further transmitter-to-receivers spacings; 

(e) means for generating a model of anisotropic formations that has horizontal resistivity Rh, vertical resistivity 
Rv, and a Dip angle with respect to a borehole reference; 

(f) means for selecting initial model values of Rh, Rv and Dip; 

(g) means for computing an error value from the differences, for each of a plurality of said transmitter-to- 
receivers spacings, between said measurement characteristics and model measurement characteristics ob- 
tained from said model using the model values of Rh, Rv, and Dip; 

(h) means for modifying said model values of Rh, Rv and Dip; 

(i) means for repeating steps (g) and (h) to minimize said error value; and 

(j) means for outputting the ultimately modified model values of Rh, Rv, and Dip. 

20. Apparatus as defined by claim 19, wherein said measurement characteristics comprise phase-shift, phase-aver- 
age, and attenuation. 

21. Apparatus as defined by claim 19, wherein said measurement characteristics comprise resistivity derived from 
phase-shift, resistivity derived from phase-average, and resistivity derived from attenuation. 

22. Apparatus as defined by claim 19, wherein said model of anisotropic formations is a transverse isotropic model. 

23. Apparatus as defined by claim 19, wherein said logging device is a logging-while-drilling device in a drill string in 
the borehole. 
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